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We present a formalism to calculate the non-linear matter power spectrum in modified gravity 
models that explain the late-time acceleration of the Universe without dark energy. Any successful 
modified gravity models should contain a mechanism to recover General Relativity (GR) on small 
scales in order to avoid the stringent constrains on deviations from GR at solar system scales. Based 
on our formalism, the quasi non-linear power spectrum in the Dvali-Gabadadze-Porratti (DGP) 
braneworld models and f{R) gravity models are derived by taking into account the mechanism 
to recover GR properly. We also extrapolate our predictions to fully non-linear scales using the 
Parametrized Post Friedmann (PPF) framework. In DGP and f{R) gravity models, the predicted 
non-linear power spectrum is shown to reproduce N-body results. We find that the mechanism to 
recover GR suppresses the difference between the modified gravity models and dark energy models 
with the same expansion history, but the difference remains large at weakly non-linear regime in 
these models. Our formalism is applicable to a wide variety of modified gravity models and it is 
ready to use once consistent models for modified gravity are developed. 

PACS numbers: 98.80.-k 



I. INTRODUCTION 



The late-time acceleration of the Universe is surely the most challenging problem in cosmology. Within the frame- 
work of general relativity (GR) , the acceleration originates from dark energy. The simplest option is the cosmological 
constant. However, in order to explain the current acceleration of the Universe, the required value of the cosmological 
constant must be incredibly small. Alternatively, there could be no dark energy, but a large distance modification 
of GR may account for the late-time acceleration of the Universe. Recently considerable efforts have been made to 
construct models for modified^ gravity as an alternative to dark energy and distinguish them from dark energy models 
by observations (see d, 1^ for reviews). 

Although fully consistent models have not been constructed yet, some indications of the nature of the modified 
gravity models have been obtained. In general, there are three regimes of gravity in modified gravity models [3, . 
On the largest scales, gravity must be modified significantly in order to explain the late time acceleration without 
introducing dark energy. On the smallest scales, the theory must approach GR because there exist stringent constraints 
on the deviation from GR at solar system scales. On intermediate scales between the cosmological horizon scales and 
the solar system scales, there can be still a deviation from GR. In fact, it is a very common feature in modified gravity 
models that there is a significant deviation from GR on large scale structure scales. This is due to the fact that, once 
we modify GR, there arises a new scalar degree of freedom in gravity. This scalar mode changes gravity even below 
the length scale where the modification of the tensor sector of gravity becomes significant, which causes the cosmic 
acceleration. 

Therefore, large scale structure of the Universe offers the best o ppo rtunity to distinguish between modified gravity 
models and dark energy models in GR 0, 0, H i, El [O, El El, 111, El, El, El El, El, [H, HH, HI . The expansion 
history of the Universe determined by the Friedman equation can be completely the same in modified gravity models 
and dark energy models. In fact, it is always possible to find a dark energy model that can mimic the expansion 
history of the Universe in a given modified gravity model by tuning the equation of state of dark energy. However, 
this degeneracy can be broken by the growth rate of structure formation. This is because the scalar degree of freedom 
in modified gravity models changes the strength of gravity on sub-horizon scales and thus changes the growth rate of 
structure formation. Thus combining the geometrical test and structure formation test, one can distinguish between 
dark energy models and modified gravity models. 
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However, there is a subtlety in testing modified gravity models using large scale structure of the Universe. In 
any successful modified gravity models, we should recover GR on small scales. Indeed, unless there is an additional 
mechanism to screen the scalar interaction which changes the growth rate of structure formation, the modification 
of gravity contradicts to the stringent constraints on the deviation from GR at solar system scales. This mechanism 
affects the non-linear clustering of dark matter. We expect that the power-spectrum of dark matter perturbations 
approaches the one in the GR dark energy model with the same expansion history of the Universe because the 
modification of gravity disappears on small scales. Then the difference between a modified gravity model and a dark 
energy model with the same expansion history becomes smaller on smaller scales. This recovery of GR has important 
implications for weak lensing measurements because the strongest signals in weak lensing measurements come from 
non-linear scales. We should note that the non-linear power spectrum is also sensitive to the properties of dark energy 

El. 

In almost all of the literature, the non-linear power spectrum in modified gravity models was derived using the 
mapping formula between the linear power spectrum and the non-linear power spectrum. This is equivalent to assume 
that gravity is modified down to small scales in the same way as in the linear regime which contradicts to the solar 
system constraints. Thus this approach overestimates the difference between modified gravity models and dark energy 
models. This was explicitly shown by N-body simulations in the context of f{R) gravity ■24. , 2 5^ , 26| . In f{R) gravity, 
the Einstein-Hilbert action is replaced by an arbitrary function of Ricci curvature (see 1271 128| for a review). This 
model is equivalent with the Brans-Dicke (BD) theory with non-trivial potential [2^l30ll3l|. The BD scalar mediates 
an additional gravitational interaction that enhances the gravitational force below the Compton wavelength of the 
BD scalar. If the mass of the BD scalar becomes larger in a dense environment like in the solar system, the Compton 
wavelength becomes short and we can recover GR. This is known as the chameleon mechanism [H. In the context 
of f{R) gravity, by tuning the function /, it is possible to make the Compton wavelength of the BD scalar short at 
solar system scales and screen the BD scalar interaction [H, [s^, [H, [s^ . N-body simulations show that, due to this 
mechanism, the deviation of the non- linear power spectrum from GR is suppressed on small scales [13, [H, [2^ . It was 
shown that the mapping formula failed to describe this recovery of GR and it overestimated the deviation from GR. 

In this paper, we develop a formalism to treat the quasi non- linear evolution of the power spectrum in modified 
gravity models by properly taking into account the mechanism to recover GR on small scales. Our formalism is based 
on the closure approximation which gives a closed set of evolution equations for the matter power spectrum [37| . 
These evolution equations reproduce the one-loop results of the standard perturbation theory (SPT) by replacing 
the quantities in the non-linear terms with linear-order ones. The SPT in GR is tested against N-body simulations 
extensively recently and it has been shown that, at the quasi non-linear regime, it can predict the power-spectrum 
with a sub-percent accuracy Although the validity regime of the perturbation theory is limited, it is the most 
relevant regime to distinguish between modified gravity models and dark energy models in GR because the difference 
in the two models is large in the linear and quasi-non-linear regime. We developed a general formalism which can be 
applied to many modified gravity models including well studied models such as f{R) gravity and braneworld models. 

This paper is organized as follows. In section II, we introduce effective equation for quasi-static perturbations 
to describe the Newtonian limit of gravity. Basic equations which are necessary to compute the power spectrum 
are presented in section III. In section IV, the non-linear evolution equations of the power spectrum are derived 
based on the closure approximation proposed by Ref . [37] . The closure approximation is one of the non-perturbative 
prescriptions for comupting non-linear power spectrum, and it is shown to be equivalent to the one-loop level of 
renormalized perturbation theory by Crocce and Scoccimarro [39] and the 2PI effective action method by Valageas 
(40| . By replacing all the quantities in non- linear terms with linear-order ones, the so-called one- loop power spectrum 
is obtained numerically which describes the leading-order non-linear corrections. In section V, we apply this formalism 
to Dvali-Gabadadze-Porratti (DGP) braneworld models. In the case of DGP models, we can derive the quasi non- 
linear power spectrum analytically. We test our numerical code to solve the closure equation against the analytical 
results. The quasi non-linear spectrum is derived also in f{R) gravity models. In this case, our results are compared 
with N-body simulations. In section VI, we apply the Parametrized post-Friedman (PPF) framework to predict the 
fully non-linear power spectrum. Using the solutions in the perturbation theory, we determine a parameter in the 
PPF framework. Then we predict the non-linear power spectrum by extrapolating this parameter to fully non-linear 
scales. The predictions of the PPF formalism are compared with N-body simulations. Section VII is devoted to 
conclusions. In appendix A, a numerical scheme to solve the closure equations is presented. In appendix B, we derive 
the quasi non-linear power spectrum in DGP models analytically. 



II. QUASI-STATIC PERTURBATIONS IN MODIFIED GRAVITY MODELS 



We consider perturbations around the Friedman-Robertson- Walker universe described in the Newtonian gauge: 

ds^ = -(1 + 2i:)dt^ + a^il + 2(j))S,jdx'dx^. (2.1) 
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We will study the evolution of matter fluctuations inside the Hubble horizon. Then we can use the quasi-static 
approximation and neglect the time derivatives of the perturbed quantities compared with the spatial derivatives. As 
mentioned in the introduction, the large distance modification of gravity, which is necessary to explain the late-time 
acceleration, generally modifies gravity even on sub-horizon scales due to the introduction of a new scalar degree of 
freedom. This modification of gravity due to the scalar mode can be described by Brans-Dicke (BD) gravity. Under 
the quasi-static approximations, perturbed modified Einstein equations give 



+ ^ = (2.2) 
2a^ 

(3 + 2cc;bd)^VV = -87rGp„(5, (2.4) 



— - 4^Gp„,(5 - — VV, (2.3) 



where G is the Newton constant measured in Cavendish-like experiments, pm is the background dark matter energy 
density and 5 is dark matter density perturbations. Under the quasi-static approximations wbd can be any function 
of time. In general, modified gravity models that explain the late time acceleration predict wbd on sub- 

horizon scales today. This would contradict to the solar system constraints which require wbd > 40000. However, 
this constraint can be applied only when the BD scalar has no potential and no self-interactions. Thus, in order to 
avoid this constraint, the BD scalar should acquire some interaction terms on small scales. In general we expect that 
the BD scalar field equation is given by 

(3 + 2cjBD)^fcV = 87rGp™(5-X(v5), (2.5) 
in a Fourier space. Here the interaction term T can be expanded as 



1 f d^ki(Pk2 
2] (27r) 
1 [ (Pki(Pk2d?k3 



= Afi(fc)(^ + - / ' Snik - fci2)Af2(fci, fc2)y^(fci)(^(fc2) 
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/ ^"^2n)f ^"" Soik - fci23)Af3(fci, fc2, k3MkiMk2Mk:,) + ... (2.6) 



where kij = ki + kj and kijk = ki + kj + fcfe. 

We should emphasize that effective equations (I2.2p . (|2.3p and (|2.5p in the BD theory can be applied only to quasi- 
static perturbations. We allow the time dependence of the BD parameter in these effective equations but this does 
not necessarily mean that we are considering scalar tensor theory where the BD parameter is a function of the BD 
scalar. The effective equations are applicable to scalar tensor theory by adding appropriate non-linear interaction 
terms in X. Although this is an interesting possibility (see [111, IH, for the analysis of perturbations in scalar tensor 
theory), we will not consider this possibility in this paper. We also note that a general parametrization of linear 
perturbations in modified gravity was developed in [4J| and a similar parametrisation of quasi-static perturbations to 
ours was considered at linearized level in [211] . 

In this paper, we consider two known mechanisms where the non-linear interaction terms X are responsible for the 
recovery of GR on small scales. There are other possibilities to recover GR on small scales for example by decoupling 
baryons, but we focus on the following two possibilities in this paper. One is the chameleon mechanism [s^ . In this 
case, the BD scalar has a non-trivial potential, and acquires a mass. Then, the BD scalar mediates the Yukawa-type 
force and the interaction decays exponentially above the length scale determined by the inverse of the mass, the 
Gompton wavelength. Because of this, the scalar interaction is hidden above the Compton wavelength, and GR is 
recovered. The BD scalar is coupled to the trace of the energy momentum tensor. Thus the effective potential depends 
on the energy density of the environment. The potential is tuned so that the mass of the BD scalar becomes large for 
a dense environment such as the solar system. In order for the chameleon mechanism to work, the scalar fields needs 
a runaway potential to be efficient ^32*1. Then the Gompton wavelength becomes very short for a dense environment 
and the scalar mode is effectively hidden. In this paper, we deal with this mechanism perturbatively. Mi determines 
the mass term in the cosmological background. The higher order terms Mi, [i > 1) describe the change of the mass 
term due to the change of the energy density. If the chameleon mechanism is at work, the effective mass becomes 
larger when the density fluctuations become non-linear. 

The other mechanism relies on the existence of the non-linear derivative interactions. A typical example is the 
Dvali-Gabadadze-Porratti (DGP) braneworld model where we are supposed to be living on a 4D brane in a 5D 
Minkowski spacetime [i^. In this model, the BD scalar is identified as the brane bending mode which describes the 
deformation of the 4D brane in the 5D bulk spacetime. The brane bending mode has a large second-order term in the 
equation of motion which cannot be neglected even when the metric perturbations remain linear. This corresponds to 
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the existence of a large M2{k) term [4y, [43, l43| ■ It has been shown that once this second order term dominates over 
the Hnear term, the scalar mode is hidden and the solutions for metric perturbations approach GR solutions. For a 
static spherically symmetric source, we can identify the length scale below which the second order interaction becomes 
important. This length scale is known as the Vainshtein radius [i^. In the cosmological situation, it is expected that 
once the density perturbations become non-linear, the second order term becomes important and we recover GR. In 
the next section, we apply the perturbation theory to solve the equations (|2.2p . (|2.3p and ()2.5|) . Thus we only keep 
up to the third order in the expansion of X which is necessary to calculate the quasi non-linear power spectrum. 

The evolution equations for matter perturbations are obtained from the conservation of energy momentum tensor, 
the continuity equation and the Euler equation; 

^ + -V-[{1 + S)v]^0, (2.7) 
ot a 

dv 1 1 

Hv + -{v ■ \7)v ^ VV". (2.8) 



at 



a a 



Eqs. (|2.3p . (|2.5p . (|2.7p and (|2.8p are the basic equations that have to be solved. In the next section, we derive evolution 
equation for perturbations in a compact form. 



III. EVOLUTION EQUATIONS FOR PERTURBATIONS 



Assuming the irrotationality of fluid quantities, the velocity field v is expressed in terms of velocity divergence 
6* = V • v/{aH). Then the Fourier transform of the fluid equations (|2.7p and (|2.8p become 



H 



dt 

_^ d0{k) 
dt 



' J "^^2^' ^°^^ ~ ~ ^(fci)<5(fe2), 



(3.1) 



H 



k 



VXfc) 



1 f dPkid^k'^ 



{27TY 



■Su{k - fci - k2)f3{ki,k2) 9{ki)e{k2), (3.2) 



where the kernels in the Fourier integrals, a and /3, are given by 



Q;(fci, fe2) = 1 



fci • k2 



/3(fcl,fc2) 



(fcl-fc2)|fcl+fc2| 
|fclP|fc2P 



(3.3) 



As for the Poisson equation (|2.3p . the potential ip is couples to S through the BD scalar in a fully non- linear way 
due to the interaction term X. To derive closed equations for 6 and 6, we must employ the perturbative approach to 
Eq. (|2.5p . By solving Eq. (|2.5p perturbatively assuming ip < 1, ip can be expressed in terms of S as 



- ^ t; 1^"^ Pr. 



3 n(fc) 



(3.4) 



where 



(3.5) 



and = SttG. The function S{k) is the non-linear source term which is obtained perturbatively using (j2.3p as 
1 f Prr,V f d^kid^k2 ^ , , ,S{ki)5{k2) 

J -wr ^'^^'^'^^''^'\{k,mk2) 



Sik) 



3 



6n(fc) 
1 

'i8n(fe) V 3 

, 6{ki)6{k2)S{k3) 



•2 - ^ f d^kid^k2d^k3 



X n I \)n.rn I I \ M2{ki,k2 + k3)M2{k2,k3) 

Snik - ki23) < Af3(A;i, A;2, fcs) — 

' ii(fes 



^23 



n(fei)n(fc2)n(fc3)' 

The expression (|3.6p is valid up to the third-order in S. 



(3.6) 
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The perturbation equations (|3.ip . (|3.2p and (|3.4p can be further reduced to a compact form by introducing the 
foUowing quantity: 



We can write down the basic equations in a single form as 

a$a(fc;T) 



m 

-9{k) 



(3.7) 



dr 



(2^) 



kid k2 ^^^^ _ ^^^^ 7af,c(fei, fc2; r) $b(fci; T)$c(fc2; r) 



(5D(fc - fci23) Cra6cd(fcl, ^2, fca; t) $6 (fcl ; t) $c (fc2T")^'d (^3 ; t), (3.. 



where the time variable r is defined by t = Ina(t). The matrix flab is given by 



^^ah(fc;r) 







V 2 



1 



1 (fc/a)2 



3 n(/s) 



-1 \ 

) 



(3.9) 



From the (2, 1) component of ilab, we can define the effective Newton constant as 

Goff = G 

If Mi=0, the effective Newton constant is given by 



l{k/al 

3 n(fc) 



2(2 + c^bd) ^ 
3 + 2cjbd 



(3.10) 



(3.11) 



See Ref. [50| for a review on the BD theory. For a positive lubd > 0, the effective gravitational constant is larger than 
GR and the gravitational force is enhanced. On the other hand, if Afi ^ fc^/a^, GeS becomes G. The quantity 'jabc 
is the vertex function as in the GR case, but new non-vanishing components arise in the case of modified gravity: 



( 1 

- a(fe2,fci) 



^abc{kl,k2;T) 



1 



a(fei, ^2) 



12i72 



'^12 



M2(fel,fc2) 



a2 ; n(fei2)n(fci)n(fc2) 



^f3{ki,k2) 







(a,6,c) = (1,1,2), 
(a,&,c) = (1,2,1), 

(a,6,c) = (2,1,1), 

(a,6,c)- (2,2,2), 
otherwise. 



(3.12) 



Here 7211 is absent in GR. Note that the symmetric properties of the vertex function, jabciki, k2', r) = 7ac6(fc2, ^i! '''), 
still hold in the modified theory of gravity. In Eq. (|3.8D . there appears another vertex function coming from the 
non-linearity of Poisson equation. The explicit form of the higher-order vertex function (Tabcd is given by 



tTa6cd(fcl,A;2,fc3;T) = < 



36i72 



K Pr, 



23 



A/3(fei,fc2,fe3 



3M3(fcl,fe2,fc3) 



n(fci23)n(fei)n(fc2)n(fc3) 

M2(fci,fc2 +fe3)M2(fc2.fc3) 



n(fc23) 



perm. 



(3 13) 

(a, 6, c, d) = (2,1,1,1'),' 
otherwise. 



Again this term is absent in GR. The vertex function (Jabcdiki, ^21 '^3! ''") defined above is invariant under the permu- 
tation of 6 <-> c <-> c? or fci <-> fc2 ^3- 
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IV. EVOLUTION EQUATIONS FOR THE POWER SPECTRUM 



In this paper, we are especially interested in the evolution of the matter power spectrum defined by 

'$a(fc;T)$fc(fc';r)\ - {2TTf S^ik + k') Pab{\k\;T). (4.1) 



Here the bracket (•) stands for the ensemble average. Note that we obtain the three different power spectra: Pss from 
(a, 5) = (1,1), -Pse from (a, 6) = (1,2) and (2,1), and Pee from (a, 6) = (2,2). 
Let us consider how to compute the power spectrum. In the standard perturbation theory (SPT), we first solve 

Eq. (|3.8|) by expanding the quantity as = ^a''' + ^'a'' + • • • . Substituting the perturbative solutions into the 
definition (j4.ip . we obtain the weakly non-linear corrections to the power spectrum. This treatment is straightforward, 
but it is not suited for numerical calculations. Furthermore, successive higher-order corrections generally converge 
poorly and SPT will be soon inapplicable at the late-time stage of the non-linear evolution. Here in order to deal 
with modified gravity models, in which analytical calculation is intractable in many cases, we take an alternative 
approach. Our approach is based on the closure approximation proposed by Ref. 37(, by which the evolution of the 
power spectrum is obtained numerically by solving a closed set of evolution equations. 

Provided the basic equation (|3.8p . evolution equations for the power spectrum can be derived by truncating an 
infinite chain of the moment equations with a help of perturvative calculations called the closure approximation. We 
skip the details of the derivation and present the final results. Readers who are interested in the derivation can refer to 
Ref. [37i] . The resultant evolution equations are the coupled equations characterized by the three statistical quantities 
including the power spectrum. We define 



$a(fc;T)$b(fc';r')) = (2^)^ fc(fc + fe') i?,,(|fe|; r, r' 
S^a{k;T 



T > T , 



<5$h(fe';r') 



= S^{k-k')Gab{k\T,T') 



T > t'. 



(4.2) 



The quantities Rab and Gab denote the cross spectra between different times and the non-linear propagator, respec- 
tively. Note that Rab 7^ Rba, in general. Then, the closure equations become 



T.abcd{k;T)Pcd{k;T) 



27r)3 



Kab{k\T) Rbcik;T, r') 



Aah{k;T) Gbc{k\T, t') 



i2n) 



(2^ 



lapq{q, k- q;T) Fbpq{-k, q,k-q:T) + -/bpqiq, k - q; r) Fapq{k, q, k - q; r) 
CTapqr (?, -q,k; T)Ppq {q; T)Prb {k;T)+ (Tbpqr {q,-q,-k;T) Ppq {q; T)Pra (fc; r) 



(4.3) 



{q, k-q]T) Kcpq{-k, q,k-q; r, r') 



dr" 
d^q 



<^apqr{q, "Q, k] T)Ppq{q; T)Rrc{k; T, t'), (4.4) 

^lapq{q,k~ q;T)jirs{-q,k;T")Gqi{\k - q\\T, T")Rpr{q;T, T")Gsc{k\T" ,t'), 



(2^)3 

d^q 
(2^3 



(q, -q, k; T)Ppg{q; T)Grc{k\T, r'). 



where the operators 'Sabcd and A^h are defined as 

d ^ d 

^abcd{k; t) = SacSbd^r + 5aAd{k] t) + Sbd^ac{k; r), Aab{k; r) = Sab^ + ^ab{k; r), 

OT OT 



(4.5) 



(4.6) 
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The explicit expressions for the kernels Fapq and Kcpq are summarized as 

Fapg{k,p,q;T) - 2 / dr" [2 G,/(g|T, r") 7irs(fc,p; r") i?ar(fc; r, r'O-Rpsfe r, r") 

+ Galik\T,T")jirsip,q;T")Rpr{p;T,T")Rg,iq;T,T")], (4.7) 

Kcpq{k',p,q;T, r') 

= 4 /"dr" G,z((7|r,r")7ir.(fe',p;r")i?p.(p;r,r") 

{Rcrik'; t', r")e(T' - r") + Rrdk'; r", r')e(r" - r')} 
+ 2 / dr"Ge/(A;'|r',r")7;rs(p,q;r")i?p.(p;r,r")i?,s(<z;r,r"). (4.8) 



To 

X 



The closure equations (|4.3P " (|4.5p are the integro-difFerential equations involving several non-linear terms, in which 
the information of the higher-order corrections in SPT is encoded. Thus, replacing all statistical quantities in these 
non-linear terms with linear-order ones, the solutions of the closure equations automatically reproduce the leading- 
order results of SPT, so called one-loop power spectra. Further, a fully non-linear treatment of the closure equations 
provides a non-perturbative description of the power spectra, and has an ability to predict the matter power spectra 
accurately beyond one-loop SPT. Strictly speaking, the non-linear terms in the right-hand side of Eqs. (|4.3|) - (|4.5[) 
have only an information of the one-loop corrections. However, it has been shown in Ref. [s?! that the present 
formulation are equivalent to the one-loop level of renormalized perturbation theory by Crocce and Scoccimarro 
(soj and 2PI effective action method by Valageas [i^, and even the leading-order approximation still contain some 
non-perturbative effects. The application of the closure approximation, together with the detailed comparison with 
N-body simulations, is presented in Refs. (ssl. [sij. Also, comprehensive discussion on the differences between several 
(semi-)analytic prescriptions for computing non-linear power spectrum, including the closure approximation, is found 
in Ref. (s^ . 

In this paper, we mainly use the closure equations for the purpose of computing the one-loop power spectra. The 
results for the fully non-linear treatment of the closure equations will be presented elsewhere. The numerical scheme 
to solve the closure equations is basically the same as the one described in Ref. In Appendix A, we briefly review 
the numerical scheme and summarize several modifications. 

Before closing this section, we note here that the resultant equations (I4.3l) - (l4.5p contain the additional non-linear 
terms originating from the modification of the Poisson equation (see Eq. (l3.4p V In particular, the terms containing 
the higher-order vertex function aabcd can be effectively absorbed into the matrix flab- Since the non- vanishing 
contribution of the higher-order vertex function only comes from a'2iii, this means that the effective Newton constant 
Gcff defined by (|3.10p is renormalised as Gcff G^s + SGcs, with SGcs given by 

3^2 r ^3g 

^^gff 47rp J (^:;;^^2iii(g,-9,fc;T)Fii(g;T). (4.9) 

This is a clear manifestation of the mechanism to recover GR on small scales through the renormalisation of the 
Newton constant and, due to this mechanism, the growth rate of structure formation is altered on non-linear scales. 



V. APPLICATIONS 



In this section, based on the formulations presented in Sec. Ulland lllll we compute the one-loop power spectrum in 
specific models of modified gravity theory. The models considered here are DGP braneworld models and f{R) gravity 
models, which we will in turn discuss in Sec. IV Al and IV Bl respectively. 



A. DGP models 

In this subsection, we consider DGP braneworld models [45j as a representative example of modified gravity models 
in the context of higher-dimensional cosmology. In this model, it is possible to derive the quasi non-linear power 
spectrum analytically as is discussed in details in Appendix B. This provides us a check of our numerical code 
explained in the previous section and Appendix A. 
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1. DGP models 



In DGP models, we are supposed to be living in a 4D brane in a 5D spacetime. The model is described by the 
action given by 

S = J d^xV^R5 + d^^V^iR + An), (5.1) 

where k"^ = SttG, R5 is the Ricci scalar in 5D and stands for the matter Lagrangian confined on a brane. The cross 
over scale Vc is the parameter in this model which is a ratio between the 5D Newton constant and the 4D Newton 
constant. The modified Friedman equation is given by 

e-=H'- ^p, (5.2) 
rc 3 

where e = ±1 represents two distinct branches of the solutions [54]. From this modified Friedman equation, we find 
that the cross-over scale Vc must be fine-tuned to be the present-day horizon scales in order to modify gravity only at 
late times. The solution with e = -f-l is known as the self-accelerating branch because even without the cosmological 
constant, the expansion of the Universe is accelerating as the Hubble parameter is constant, H ~ l/vc- On the other 
hand e = — 1 corresponds to the normal branch. In this branch, we need a cosmological constant to realize the cosmic 
acceleration. However, due to the modified gravity effects, the Universe behaves as if it were filled with the Phantom 
dark energy with the equation of state w smaller than —1 [5^ [5^. It is known that the self- accelerating solution is 
plagued by the ghost instabilities (see f57l| for a review). Also it gives a poor fit to the observations such as supcrnovae 
and cosmic microwave background anisotropics [58, 59 , 60.] . 

However, this model is the simplest modified gravity model where the mechanism to recover of GR on small scales 
is naturally encoded and it can be used to get insights into the effect of this mechanism on the non-linear power 
spectrum. Also in this model, it is possible to derive analytic expressions for the quasi non-linear power spectrum. 
Thus it offers a check of our numerical code to solve the evolution equations for the power spectrum. 

In this model, gravity becomes 5D on large scales larger than rc. On small scales, gravity becomes 4D but it is not 
described by GR. The quasi-static perturbations are described by the BD theory where the BD parameter is given 
by ^ 

^BD(T) = ^(/3(r)-l), (3{T) = l-2eHrc(^l + ^y (5.3) 

where H is the cosmic time derivative of the Hubble parameter H . Note that the BD parameter depends on time in 
this model. The BD scalar is massless A/i = 0. However, it acquires a large second order interaction given by [4^ 



I(^) 



(V2^)2 - (V,Vj^)2 . (5.4) 



Note that Tc is tuned to be the present-day horizon scale. Thus this second order term has a large effect. The 
higher order terms than the second order are suppressed by the additional powers of the 4D Planck scale and, in the 
Newtonian limit, we can safely ignore them. Therefore, in this model, we have 

Mi=0, M2{ki,k2) ^ 2^ kjk^-- (ki-k2)\ M3 = 0, (5.5) 

Uik,T) = /3(t)^. (5.6) 

In this paper, we use the best-fit cosmological parameters for the flat self-accelerating universe: flm — 0.258, ilb = 
0.0544, h = 0.66, — 0.998 [60]. For the normal branch, we add the cosmological constant fl^ = 1.5. 

It has been suggested that the idea of the self-acceleration could be extended to more general models and in these 
models one could avoid the ghost instabilities [62] . Although the covariant theory that realizes this idea is not known 
yet but our formalism can be applied to calculate the non-linear power spectrum in these models. In fact, in our 
formalism, this new model corresponds to add a new term M3 (and also a constant term in X). We will leave it as a 
future work to study these extensions. 
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FIG. 1: Fractional change in the non-linear power spectrum relative to the linear power spectrum in the self-accelerating branch 
of DGP. The solid (black) line is the analytic solutions and the circles are numerical solutions obtained by solving the closure 
equation. The dashed (red) line shows the analytic solutions obtained by neglecting the non-linear interaction terms I. The 
triangles represent the numerical solutions in this case. We used the best fit cosmological parameters for the flat universe 
an = 0.258, Qb = 0.0544, h = 0.66, Us = 0.998. 



Fig. 1 sliows a fractional cliange in the non-linear power spectrum for density perturbations relative to the linear 
power spectrum in the self-accelerating branch. In order to see the effect of the non-linear interaction term X, we also 
plotted the non-linear power spectrum obtained from the linear Poisson equation by neglecting the interaction term 
X. Although the differences between the two spectra are small in the quasi non-linear regime, we can see that the 
non-linear interaction term X enhances the non-linear power spectrum. This is natural because in the self-accelerating 
branch, the linear growth rate is suppressed compared to the GR model that follows the same expansion history. The 
suppression is due to the negative BD parameter lobb < in this model, which makes the Newton constant smaller 
than GR. This is closely related to the fact that the BD scalar becomes a ghost. Classically, the ghost mediates a 
repulsive force and suppresses the gravitational collapse. The non-linear interaction makes the theory approach GR. 
Thus it effectively increases the Newton constant by screening the BD scalar. Then the power spectrum receives an 
enhancement compared with the case without the non-linear interaction. 

Fig. 2 shows a fractional change in the non-linear power spectrum relative to the linear power spectrum in the 
normal branch. Again we showed the two cases with and without 2. In the normal branch, the linear growth rate is 
enhanced as the BD parameter is positive wbd > 0. Thus, the situation is completely opposite to the self- accelerating 
branch and the non-linear interaction suppresses the non-linear power spectrum in order to recover GR on small 
scales. 

In both branches, it is possible to derive the solutions for the power spectrum analytically by solving the equation 
for perturbations (j3.8p perturbatively 



2. Quasi non-linear power spectra 



(5.7) 



The power spectrum is also expanded accordingly 



(5.8) 



The detailed calculations are summarized in Appendix B. In Fig.l, we compare the results obtained by solving 
the closure equations numerically with those from the analytic solutions. In order to derive the analytic solutions. 



10 




k \h Mpc-^] 



FIG. 2: The same in the normal branch as in Fig.l. We only show the analytic solutions. The numerical solutions agree 
with them very well. The cosmological parameters are the same as the self-accelerating universe but in addition there is a 
cosmological constant SIa — 1.5. 



we employed the Einstcin-de Sitter (EdS) approximation. In the EdS approximation, aU the non-linear growth 
rates appearing in the higher-order solutions are approximately determined by the linear growth rate Di(t). It is 
also possible to apply the EdS approximation in the numerical calculations [53| and we have checked that the EdS 
approximation changes the result only at sub percent level. The fact that the two results agree very well confirms the 
validity of our numerical code. 



B. f{R) gravity models 

In this subsection, we derive the quasi non- linear power spectrum in f{R) gravity model (see (27l. [28j for reviews). 
In this model, N-body simulations have been performed [2^l25l [2^ and we will check our numerical solutions against 
the full N-body simulations. 



1. f{R) gravity models 

We consider another class of modified theory of gravity that generalizes the Einstein-Hilbert action to include an 
arbitrary function of the scalar curvature R: 



R + fjB) 
2k2 



in 



(5.9) 



where — SttG and L^i is the Lagrangian of the ordinary matter. This theory is equivalent to the BD theory with 
i^BD = but there is a non-trivial potential (63l. [6^. This can be seen from the trace of modified Einstein equations: 



(5.10) 



where fn = df /dR and □ is a Laplacian operator and we assumed matter dominated universe. We can identify fji as 
the BD scalar field and its perturbations are defined as 



(5.11) 
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where the bar indicates that the quantity is evaluated on the cosmological background. In this paper, we assume 
I/rI ^ 1 and \ f/R\ ^ 1- These conditions are necessary to have the background close to ACDM cosmology. Then 
the BD scalar perturbations satisfy 

S^VV^-'^W + 'Ji?, 5R=R{fn)~RCfR). (5.12) 

This is nothing but the equation for the BD scalar perturbations with t^BD = and the potential gives the non-linear 
interaction term 



Then we find 



J(^) = 5R{lp). (5.13) 



kV , Rf{r) 



n(fc,T)= - 



a 

We should note that in this model, the linear growth rate depends on the wave number. Due to this, the vertex 
functions are not the separable functions in terms of k and r. This prevents us deriving the solutions analytically 
unlike the DGP case and we need to solve the closure equation directly. 
In this paper, we consider the function f{R) of the form 

f{R)(x ^ ^ , (5.14) 

where ^4 is a constant with dimensions of length squared [s^. In the limit R ^ 0, f{R) and there is no 
cosmological constant. For high curvature AR 3> 1, f{R) can be expanded as 

/(i?)~-2KVA-/™^, (5.15) 

where pA is determined by A, Rq is the background curvature today and we defined ffiQ as //jo — fniRo)- As we 
mentioned before, we take |/fl,o| ^ 1 and assume that the background expansion follows the A CDM history with the 
same pA- The Mi term determines the mass of the BD field ttibd — {Mi/Sy^^ as 




(5.16) 

Above the Compton length nij^^, the BD scalar interaction decays exponentially and we recover GR. On small scales, 
we recover the BD theory with cjbd = if wc neglect the higher order terms Ali, {i > 1). From Eq. (13. lip , the Newton 
constant is 4/3 times large than GR. Thus the linear power spectrum acquires a scale dependent enhancement on 
small scales. Of course, this model is excluded from local gravity constraints. The higher order terms Mi{i > 1) are 
responsible for suppressing this modification of gravity on small scales via the chameleon mechanism and it makes 
possible to pass local gravity constraints. Thus we expect that the non-linear interaction terms 2 will suppress the 
non-linear power spectrum. In the following, when we mention the power spectrum with the chameleon mechanism, 
it means that we introduce the non-linear terms X derived from (|5.14p which contains the mechanism to recover 
GR on small scales by screening the scalar mode. In this paper, we adopt the cosmological parameters given by 
l/flol = 10-4, ris = 0.958, = 0.24, = 0.046, Qa = 0.76, h = 0.73. 

In this model, the solar system constraints are satisfied but it has be en p ointed out that the chameleon mechanism 
does not work for strong gravity and neutron stars cannot exist [1^ |66| |. It was claimed that a fine-tuned higher 
curvature corrections to f{R) is needed to cure the problem fEf]. However, recently there appeared an objection 
against the absence of relativistic stars [68| . In this paper, we perturbatively take into account the chameleon 
mechanism in the cosmological background and the quasi non-linear power spectrum would be insensitive to the 
higher curvature corrections. 

We should also mention that recently a couple of papers appear that study the second order perturbations and 
three point functions in f{R) gravity models [69l. F70l|. 
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FIG. 3: Fractional change in the non-hnear power spectrum in f{R) gravity models relative to the GR models with the same 
expansion history. Pgr is the non-linear power spectrum in ACDM model with the same cosmological parameters. The solid 
(black) lines are the solutions in the perturbation theory obtained by solving the closure equation numerically. The circles 
show the results of N-body simulations. The dashed (red) line is the perturbation theory solutions obtained by neglecting the 
non-linear interaction terms T. The triangles represent the corresponding N-body solutions. The arrow indicate the valid regime 
of the perturbation theory. The parameters are taken as \ fRo\ = 10"*, n, = 0.958, fl^ = 0.24, fib = 0.046, = 0.76, h = 0.73. 



2. Quasi non-linear power spectra 



Fig. 3 shows a fractional change in the non-hnear power spectrum in f{R) gravity models relative to the GR model 
with the same expansion history, ACDM model. Note again that unlike the DGP case, the function 11 is not a 
separable function of k and t, and the analytical calculation is intractable. The solid line shows the case with the 
chameleon mechanism which includes the higher order terms M2 and M3. The dashed line shows the non-chameleon 
case where we neglected the higher order terms A/2 and M3. As expected, the chameleon mechanism suppresses the 
non-linear power spectrum. The situation is similar to the normal branch DGP models. 

We have also shown the results from N-body simulations in both cases [7l| . The N-body simulations are performed 
under the quasi-static approximations. The form of f{R) that is used in the simulations is the same as ours and we 
have checked that the basic equations that are used in their simulations are the same as our equations. They run 
the simulations with a box size 256h~^ Mpc and with 256^ particles. The power spectrum starts to show systematic 
> 10% deviations from Smith et.al. fitting formula [tI for fc > 0.79/i Mpc The arrow in the figure shows the 
expected validity range of the perturbation theory. We determine this validity regime of the perturbation theory by 
solving the equation 

f d^qPUq,z) ^0.18, (5.17) 
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for k, where Pun is the hnear power spectrum. This Hmit of k represents the range of the 1%-level accuracy, which 
has been empirically found by comparisons between the perturbation theory predictions and N-body simulations [38j . 
Of course, this condition was calibrated in simulations in GR and there is no guarantee that this condition can be 
applied to modified gravity models, but this limit gives an useful indication of the validity of the perturbations theory. 
In fact, we find that the agreement between the perturbations theory and N-body simulations in f{R) theory is good 
within this validity regime of the perturbation theory. 



VI. IMPLICATION FOR FULLY NON-LINEAR POWER SPECTRUM 



We have developed the formalism to derive the quasi non-linear power spectrum. In this section, we study the 
Parametrized Post-Friedmann (PPF) framework for the non-linear power spectrum Q in order to get insight into the 
fully non-linear power spectrum. 



A. PPF formalism 

Hu and Sawicki proposed a fitting formula for the non-linear power spectrum in modified gravity models [s']. The 
fitting formula based on the observation that the non-linear power spectrum should approach the one in the GR model 
that follows the same expansion history of the Universe due to the recovery of GR on small scales. They postulate 
that the fully non-linear power spectrum in a modified gravity model is given by the formula 

p/, ^ _ Pnon-GR{k,z)+Cnl^'^{k,z)PQnik,z) 

l + c,,E2(fc,z) ' ^^-^^ 

where z is a red-shift. Here P^ion-GRik, z) is the non- linear power spectrum which is obtained without the non-linear 
interactions that are responsible for the recovery of GR. This is equivalent to assume that gravity is modified down 
to the small scales in the same way as in the linear regime. Pcnik, z) is the non-linear power spectrum obtained in 
the GR dark energy model that follows the same expansion history of the Universe as the modified gravity model. 
The function Y?{k,z) determines the degree of non-linearity at a relevant wavenumber k. They propose to take 
I]^(fc) — k^P\in{k, z)/2tt'^, where Piin(^, z) is the linear power spectrum in the modified gravity model. Finally, Cni is 
a parameter in this framework which controls the scale at which the theory approaches GR. 

Once we obtain the quasi non-linear power spectrum, we can check whether the PPF framework works and deter- 
mines Cni in the quasi non-linear regime. In our formalism, Pnon-GR{k, z) is obtained by neglecting the non-linear 
interaction X. PGR(fc, z) can be obtained by taking ujbd oo limit and also neglecting X. We again consider two 
explicit examples. 



B. DGP models 



We first consider the self-accelerating branch solutions in DGP models. The extension to the normal branch 
solutions is straightforward though there is a subtlety in defining the equivalent dark energy model as the equation 
of state of dark energy becomes less than —1 and it diverges at some redshift [t^. Also an addition of curvature 
in the background is necessary to have large modified gravity effects [74|. We leave this for a future work. In the 
self-accelerating universe, we find that 

S2(fc,z) = ^Pun(fc,z), (6.2) 

as proposed by Hu and Sawicki Q gives a nice fit to the results obtained in the perturbation theory. We find that 
by allowing the time dependence in Cni, it is possible to recover the solutions for the non-linear power spectrum very 
well within the validity regime of the perturbation theory determined by Eq. (|5.17p . At z = 0, Cni is given by 0.3 
and there is a slight redshift dependence (Fig. 4). For < z < 1, we find that Cni can be approximately fitted as 
c„i = 0.3(1 -h 

Armed with this result, it is tempting to extend our analysis to the fully non-linear regime. In GR, there are several 
fitting formulae which provide the mapping between the linear power spectrum and non-linear power spectrum. It is 
impossible to apply these mapping formulae to modified gravity models as the mapping does not take into account 
the non-linear interaction terms X in the Poisson equation. If we apply the GR mapping formula to the linear power 
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spectrum, we would get the non-linear power spectrum without X, that is, Pnon-GR(fc)- In fact, there exists N-body 
simulations in DGP models performed by neglecting the non-linear interaction terms. 

The mapping formulae should be valid in GR models, so we can also predict PgrC^)- Then using the PPF formalism 
(|6.ip . we can predict the non-linear power spectrum if c^i is known. In the right panel of Fig. 4, we plotted a fractional 
change in the power spectrum in the DGP model relative to the GR model with the same expansion history. We used 
the fitting formula developed by Smith et.al. It should be noted that this fitting formula should be used with 

caution even in the simple case of dynamical dark energy 75] but given a lack of better fitting formula available at 
the moment, we used this formula as an example. If we could extrapolate the result in the quasi non-linear regime, 
we would have Cni = 0.3 at z = 0. 

Recently, N-body simulations in DGP models have been done by Schmidt JTg]. Fig. 5 shows the comparison between 
the PPF prediction and N-body simulations. In the left panel of Fig. 5, the dashed line shows the power spectrum with 
the linear Poisson equation without X. The corresponding N-body results are shown by triangles. Although Rcf. [T^ 
reported that the fitting formula works fine in this case, we find that the fitting formula by Smith. et.al. slightly 
overestimates the power spectrum [tHI • The solid curve shows the full power spectrum including 1. Again the PPF 
formalism slightly overestimates the power spectrum. If we take the ratio between the two cases, the PPF formalism 
recovers N-body results up to fc = O.SMpc h~^. The validity regime of perturbation theory is below k = 0.12Mpc 
h"^. Thus the PPF formalism using Cni derived by the perturbation theory describes the effect of the Vainshtein 
mechanism on non- linear scales beyond the validity regime of the perturbation theory. 

This observation suggests that an improvement of the PPF formalism can be made by getting a more accurate 
power spectrum without the Vainshtein mechanism because the PPF formalism describes the effect of the Vainshtein 
mechanism very well. In order to demonstrate this fact, we derive the power spectrum without the Vainsgtein 
mechanism Pnon-GR by interpolating the N-body results (see the right panel of Fig. 5). Using this power spectrum 
as the power spectrum with the linear Poisson equation Pnon-GR in the PPF formalism, we find that the full power 
spectrum can be very well described by the PPF formalism where Cni is derived by the perturbation theory. We 
should emphasize that the ratio between the power spectra with and without the chameleon mechanism is not very 
sensitive to the power spectrum with the linear Poisson equation. This also indicates that the PPF formalism with 
Cni determined by the perturbation theory describes the effect of the Vainshtein mechanism very well at least up to 
k - 0.5/iMpc"^ 

Above k = IMpch^^, N-body simulations do not have enough resolutions. If we can extrapolate our results to 
this regime, we find that even aX k = lOMpch^^, the difference between the power spectrum in DGP and that in the 
equivalent GR model remains at 7% level. This is crucial to distinguish between the two models using weak lensing 
as the signal to noise ratio is larger on smaller scales. Of course, we should emphasize that there is no guarantee that 
Cni measured in the quasi non-linear regime is valid down to the fully non-linear scales and this should be tested by 
N-body simulations with higher resolutions. 



C. f(R) gravity models 

Next, we consider f{R) gravity models. As in the DGP model, we first check if we can reproduce the perturbation 
theory results by the PPF fitting. We find that the fitting is not very good if we adopt S(fc, z) as is proposed by Hu 
and Sawicki [1]. Instead, if we choose I](fc, z) as 

/ fc3 \ 1/3 

^^ik,z)^ i^PUk,z)\ , (6.3) 

the solutions in the perturbation theory are fitted by the PPF formalism very well by allowing the redshift dependence 
in Cni. At z = 0, Cni = 0.085 gives an excellent fit to the power spectrum within the validity regime of the perturbation 
theory. For < z < 1, we find that Cni can be approximately fitted as Cni = 0.08(1 -I- z)! "^. In Fig. 6, we also show 
the prediction for the fractional difference between the power spectrum in f{R) theory and that in ACDM model in 
fully non- linear regime for several Cni. 

It is also possible to check our predictions against the N-body simulations done by Refs. [13, H^- Fig. 7 
shows the comparison between the PPF prediction and N-body simulations. In the left panel, the dashed line 
corresponds to non-chameleon case with Cni = 0. The corresponding N-body results are shown by triangles. We again 
used the fitting formula by Smith et.al. to derive the non-linear power spectrum from the linear power spectrum. 
Compared with the N-body results, we find that, in this case, the formula by Smith et.al. slightly underestimates 
the power spectrum around 0.03/iMpc"i < k < 0.5hMpc~^ and overestimates the power at fc > 0.5hMpc~^ though 
N-body simulations have large errors in this regime. The solid line shows the case with the chameleon mechanism. 
Again the PPF formalism underestimates the power spectrum in the same region as the non-chameleon case. If we 
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FIG. 4: A fractional change in the power spectrum in the DGP self-accelerating solution relative to the GR model which has 
the same expansion history as the DGP. The solid (black) line shows the perturbation theory solution and the dashed (red) line 
shows the perturbation theory solution without the non-linear interaction terms in the Poisson equation. The dotted (blue) 
line shows the PPF fitting. By allowing the redshift dependence of Cni, we can fit the power spectrum very well within the 
validity regime of the perturbation theory indicated by arrows. The right panel shows the results at 2 = obtained from the 
fitting formula by Smith et.al. for Pnon-GR and Pgr- If Cni ~ 0.3 obtained by the perturbation theory is applicable, the solid 
(black) line is our prediction on non-linear scales. The cosmological parameters are the same as in Fig 1. 



take the ratio between the non-chameleon case and chameleon case, the PPF formalism nicely recovers the N-body 
results up to fc ~ 0.5ft.Mpc^^. Beyond that, N-body simulations have large errors. We should emphasize that the 
perturbation theory is valid only up to fc = 0.08/iMpc^^ at z = 0. Thus the PPF formalism using c^i derived by the 
perturbation theory describes the effect of the chameleon mechanism on non-linear scales beyond the validity regime 
of the perturbation theory. 

As we have done in DGP models, we also derived the power spectrum without the chameleon mechanism Pnon-GR 
by interpolating the N-body results (see the right panel of Fig. 7). Using this power spectrum as the non-chameleon 
power spectrum Pnon-GR in tke PPF formalism, we find that the power spectrum with the chameleon mechanism 
can be very well described by the PPF formalism where Cni is derived by the perturbation theory. Again for larger fc, 
N-body simulations also do not have enough resolutions and it is difficult to tell whether this extrapolation is good 
or not. More detailed study is needed to address the power spectrum at larger k, but the PPF formalism is likely to 
give a promising way to develop a fitting formula for the non-linear power spectrum in modified gravity models. 



VII. CONCLUSIONS 



In this paper, we developed a formalism to derive the quasi-nonlinear power spectrum for density perturbations 
in modified gravity models. We assume that the quasi-static perturbations under the horizon scale can be described 
by a BD theory with a time dependent BD parameter. The non-linear interaction terms are model-dependent so 
we parametrised these terms. Then we took into account the non-linear interaction terms perturbatively. The 
Poisson equation becomes a non-linear equation which relates the curvature perturbation to the density perturbations 
non-linearly. Then combining it to the continuity equation and Euler equation assuming the irrotationality of fluid 
quantities, we derived the evolution equations for the density perturbations and the velocity divergence in a compact 
form. The non-linearity of the Poisson equation introduces new vertex functions. 

Then we derived the evolution equations for the power spectrum. The closure approximation was employed to 
derive a closed set of equations. In this paper, we further simplified the equations by replacing all the quantities 
in non-linear terms with linear-order ones. The resultant theory is equivalent to the standard perturbation theory 
where we solve the perturbations up to the third order. The advantage of using the closure equations is that we can 
directly integrate the evolution equations for the power spectrum numerically without using approximations to derive 



16 



-0.1 
-0.15 

-0.2 

-0.25 

-0.3 
1.1 

1.05 




A linear Poisson eq. t j 
• sDGP 



Smith et.al. 

Lines: PPF 
Symbols: N-body 



H 1 1 I I I I I 




-0.1 



fi -0.15 



o 



O 



-0.2 



-0.25 



-0.3 



1.1 



e 1 

0.01 



2=0 sDGP 








\ 

\ 

X 


- 


A linear Poisson eq. 


St 

V : 


• sDGP 


fitting to linear Poisson eq. 


\ 


Lines: PPF 


\ 


Symbols: N-body 


\ 


1 1 1 1 — 1 — 1 1 1 1 1 1 h- 


-y- 

a/ 



0.1 1 

k [h Mpc^] 

FIG. 5: Comparison between the PPF prediction and N-body simulations. In the left panel, Smith et.al. fitting formula is 
used to predict Pnon-GR and Pgr- We used Cni determined by the perturbation theory Cni = 0.3 at z = 0. In the right panel, 
we fitted N-body results with the linear Poisson equation to derive Pnon-GR- 



solutions for perturbations. The analysis of the full closure equations will be presented in a separate publication. 

We solved the closure equations in DGP and f{R) gravity models as examples. In the DGP model, the two 
branches of the solutions were considered. In the self-accelerating branch, the expansion of the Universe is accelerated 
without having the cosmological constant. In this branch, the BD parameter is negative and the linear growth rate 
is suppressed compared with the dark energy model that follows the same expansion history as the self-accelerating 
universe. We found that the non-linear interactions in the Poisson equation recover GR hence enhance the power 
spectrum by shielding the BD scalar interactions. On the other hand, in the normal branch solutions, where we need 
the cosmological constant to accelerate the expansion of the universe, the BD parameter is positive. Then the linear 
growth rate is enhanced compared to the corresponding dark energy model. In this case, the non-linearity of the 
Poisson equation suppresses the power spectrum. In DGP models, it is also possible to derive the semi-analytical 
solutions for the power spectrum by using the Einstein-de Sitter approximations for the growth rate. It was shown 
that the two methods give almost identical results. In /(i?) gravity models, it is impossible to derive semi-analytic 
solutions as the linear growth rate depends on scales. In this case, we compare our results with N-body simulations. 
Within the validity regime of the perturbation theory inferred from an empirical formula derived in GR simulations, 
our numerical solutions agree with N-body results very well. In f{R) gravity models, the BD parameter is vanishing 
and gravity becomes strong below the Compton wavelength of the BD field and the linear growth rate acquires a 
scale dependent enhancement. Then the chameleon mechanism that recovers GR on small scales suppresses the power 
spectrum. 

Recently, the parametrization of the non-linear power spectrum in modified gravity models was proposed within 
the framework of Parametrized Post Friedman (PPF) formalism. We checked whether the PPF formalism works or 
not explicitly using the solutions in the perturbation theory. In DGP models, we find that the PPF formalism works 
very well within the validity regime of the perturbation theory by allowing a time dependence in the PPF parameter 
Cni. The power spectrum without the mechanism to recover GR on small scales can be predicted by using a mapping 
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formula between the linear power spectrum and non-linear power spectrum derived in GR. We used a fitting formula 
obtained by Smith et.al which gives an accurate non- linear power spectrum in GR. Using the value of Cni obtained in 
the perturbation theory, we predicted a fractional change in the non-linear power spectrum in DGP relative to the 
equivalent dark energy models. With the best fit cosmological parameters for SNe and CMB, there is a 14% change 
on linear scales and there is still a 7% change at fc = 10/iMpc~^ at z = 0. Interestingly, the maximum change appears 
around k = l/iMpc~^. We apphed the same procedure to f{R) gravity models. It was found that a slight modification 
was needed in the definition of z) in the PPF formalism to reproduce the solutions in the perturbation theory. 
With this modification, we could reproduce the perturbation theory solutions very well again by allowing a time 
dependence in c^i- Then we predicted the fully non-linear power spectrum again using Smith et.al. fitting formula 
and the value of Cni obtained in the perturbation theory. 

We compared our PPF predictions with the N-body results in both cases. It was found that the PPF formalism 
reproduces the ratio between the power spectrum with and without the Vainshtein/chamelcon mechanism very well 
beyond the validity regime of the perturbation theory. This indicates that the PPF formalism calibrated by the pertur- 
bation theory describes the effect of the Vainshtein/chameleon mechanism very well. On the other hand. Smith et.al. 
fitting formula does not predict the power spectrum without the Vainshtein/chameleon mechanism very accurately. In 
f{R) gravity, this would be because the linear growth rate depends on scales. An improvement of the prediction can 
be made by obtaining the power spectrum without the Vainshtein/chameleon mechanism more accurately. In order 
to demonstrate this fact, we fitted the N-body results without the Vainshtein/chameleon mechanism and applied the 
PPF formalism to predict the power spectrum with the Vainsgtein/chameleon mechanism. In this way, it was shown 
that the PPF formalism calibrated by the perturbation theory reproduces the N-body results within the errors in 
N-body simulations. It should be pointed out that N-body simulations without the Veinshtein/chameleon mechanism 
is much easier to perform as the modified gravity effects reduce to the change of the Newton constant. On the other 
hand, if we need to include the Veinshtein/chameleon mechanism properly, it is necessary to solve the Klein-Gordon 
equation for the scalar field directly. 

The understanding of the non-linear power spectrum is essential to distinguish between modified gravity models 
and dark energy models. Especially, weak lensing is sensitive to the clustering property on non-linear scales and all 
the predictions done by neglecting the mechanism to recover GR on small scales should be revisited. Our formalism 
based on the perturbation theory enables us to predict the quasi non-linear power spectrum very accurately. The 
PPF formalism calibrated by the perturbation theory gives an analytical way to predict non-linear power spectrum. 
Of course, we eventually need N-body simulations to check the predictions, but the PPF formalism will provide a way 
to develop a fitting formula for the non-linear power spectrum which is widely used to predict weak lensing signal in 
GR dark energy models. We tested our formalism in DGP and f{R) gravity models but we should bear in mind that 
these models face serious difficulties. Our formalism is applicable to a wide variety of modified gravity models and it 
is ready to use once consistent models for modified gravity are developed. 
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FIG. 7: Comparison between the PPF prediction and N-body simulations in f{R) gravity models. In the left panel, Smith 
et.al. fitting formula is used to predict Pnon-GR and Pgr. We used Cni determined by the perturbation theory Cni — 0.085 at 
z = 0. In the right panel, we fitted N-body results without the chameleon mechanism to derive Pnon-GR. 
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APPENDIX A: NUMERICAL ALGORITHM 

In this paper, we follow the numerical scheme presented in [53| to solve the closure equations derived in [33]. The 
basic procedure to solve Eqs. (|4.3p - (|4.5p is the same as the one summarised in Sec. Ill in [ssj. In what follows, we 
focus on the modifications that are necessary to apply it to modified gravity models. 

In (53| . the integration over k has been performed before performing the time evolution thanks to the time- 
independent vertex functions which make the integrands separable functions in terms of time and k. In the present 
case, however, the integration over k has to be performed at each time step because the integrands are generally not 
separable due to the time-dependent vertex functions. Another different point from the previous work is the existence 
of the extra vertex function, cr2iii. 

Here we extend our numerical scheme presented in [s^ to include the above two new ingredients, the non-separable 
integrands and the extra vertex function. Firstly, we change the variable k' to k — k' in the integrations over k' in 
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Eqs. (|4.3p - (|4.5p . then we obtain a more symmetric form of the closure equations, 

KbGbMr, t') = r dr" Masik; r, T")G,,(fc|T", r') + Sar{k; r)G,,(fc|T, r'), (Al) 



where 



AatRUk;T,T')= / dr" Mas{k;T,T")Rsc{k;T",T')+ dr" Nai{k;T,T")G,e{k\T' ,t") 

J To J To 

+ Sar{k;T)Rrc{k;T,T'), (A2) 

tabcdPcd{k;T) ^ / dr" Mas{k;T,T")Rbs{k;T,T") + / dr" Nae{k;T,T")Gb£{k\T,T") 

J To J To 

+ Sar{k;T)Prbik;T) + {a^b), (A3) 

Mas{k;T,T")=A j -^^^p^{k- k\k';T)-fersik' - k,k;T")G,eik'\T,T")Rpri\k- k'\;T,T"), (A4) 

Naiik;r,T")^2 j ^^7ap,(fc - fc', fc'; r)7fr,(fe - fe', fe'; r")i?,.(fc'; r, T")^pr(|fc - fe'l; r, r"), (A5) 
d^k ' 

W) 



Sar{k;T) = 3 / -j—-:^aapqr{k',-k',k;T)Ppq{k';T), (A6) 



and 



The domain of the two-dimensional integrand (|A4p - (jA6p is determined by 

^min !i 1^ I: 1^ k \ < fcmax: (A8) 

where we set fcmin = O.OOl/iMpc^^ and fcmax — 5/iMpc~^ according to the convergence check performed in [5^. 

Next we introduce {X, Y, (/))-coordinates in the k' space. Suppose k is aligned to the k'^ axis in fc' space. Then the 
elliptic coordinate is defined as 



^ / sinh sin /z cos ( 



k' = — + q, Q=\ sinh C sin /i sin . (A9) 
\ cosh ( cos fi 



X^coshC, r = cos/x, (AlO) 



2 

We introduce the XF-coordinate as 



where X > 1 and — 1 < < 1. In these coordinates, the domain (|A8|) is depicted as the shaded hexagonal shape in 
FigEl Then we perform the integration over the domain at each time step. Note that the integration over yields only 
a constant 27r. In the previous work, we used the trapezium rule for the integrations. Instead, in order to reduce the 
computational cost without loosing the numerical accuracy, we implement the two-dimensional Gaussian quadrature 
(GQ). The two-dimensional GQ is adopted to the circumscribed rectangle of the computational domain as shown in 
Fig. ([8]). The circumscribed rectangle is then transformed to a unit square domain, {x,y) G [— 1 : 1] x [— 1 : 1]. On 
the unit domain, the two-dimensional integration is approximated by 

»1 »i Ngq 

/ / fix,y)dxdy K f{xi,y^)wiWj, (All) 
J-iJ-i 

where Xi and yi are defined as zero points of the Legendre function, namely, Pngq i^i) = 0. The weight Wi is calculated 
by 

(A12) 
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where the prime denotes derivative with respect to x (e.g., see [78|). The domain of the two-dimensional integrands 
(|A8[) is extended to the unit domain so that the integrands vanish on the extended region, and those on the grid 
{xi , yj ) are computed by the cubic sphne interpolation of the propagators and power spectra. Having checked the 
convergence of the numerical results, we fix the number of the grid points as Nqq — 31. 

For the time evolutions, we set the initial redshift as zq = 400, and the number of time steps as ~ 400, which 
are chosen so that the final results become insensitive to these parameters at the order of 0.1%. The number of wave 
numbers for the propagator and power spectra is set to be = 100, and the discrete points in the wave number 
space are taken so that they become dense around the baryon acoustic oscillation scale. 




FIG. 8: The domain of the integrations in the integrations (|A4[) - (|A6[) in the (X, y)-coordinate. The shadow region is defined 
by Eq. (|A8|I and, by definition, — 1 < ^ < 1. The circumscribed rectangle is the computational domain for the Gaussian 
quadrature. 



APPENDIX B: PERTURBATION THEORY AND ONE-LOOP POWER SPECTRA IN DGP MODEL 



In this appendix, we derive a set of perturbative solutions up to the third order in DGP models. Based on these 
solutions, we obtain the analytic expressions for the one- loop power spectra. 



1. Solutions up to the third order 



The solutions for 5 and 9 are perturbatively expanded as 

5{k) = (5(i)(fe) + S'^^Kk) + (5(3^^) + • • • , S{k) ^ 6'(i)(fc) + + 6'(3)(fc) + . . . . (Bl) 

In order to obtain the solutions analytically, we use the Einstein-de Sitter (EdS) approximation (see sec. 2. 4. 4 in [t^ 
for a review). In the EdS approximation, all the non-linear growth rates appearing in the higher-order solutions are 
approximately determined by the linear growth rate Di{t). The explicit expressions for solutions at each order can 
be obtained analytically, and are summarized below. 



Ist-order solutions : 

<5W(fc;t) = Di(t)<5o(fc), B^^\k;t) = -{D^/ H)5^{k), (B2) 

with Di being the linear growth rate. Here a dot denotes the derivative with respect to time and (5o(A;) denotes the 
primordial density perturbation defined at early matter dominated era. We assume 5{){k) obeys Gaussian statistic. 
The evolution equation for Di is given by 



CDx = 0, 



(B3) 
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where the linear operator C, is given by 



^ d2 d K ( \ 



2nd-order solutions: 



Su{k - fci2) (5o(fci)^o(fc2) [Dlit) F^p^{kl,k2) + F2{t){l - fil^) 

r rl^k rl^kr, r 
0^^\k;t) = I ^ fe(fc - fci2) ^o(fci)^o(fc2) -(Di£'i//f)GgL(fei,fe2)-(i^2/if)(l-Mf,2) 



(27r)3 



where we define 



Mi 



\k,\\k,\- 



(B4) 

(B5) 
(B6) 

(B7) 



The symmetrized kernels of the integrals, -Fgym and Gsym, are respectively defined as 

F^^l{ki,k2) = ^{a(fei,fe2) + a(fe2,fei)} + ^/3(fei,fe2), 

3 2 
Gg|„(fei,fe2) = Y^{aiki,k2) + a{k2,ki)} + -/3{ki,k2). 

Note that in addition to the hnear growth rate Di, the second-order solutions contain the new growth function F2, 
which originates from the non-hnearity of the Poisson equation. The evolution equation for F2 is 



rP — ( I^Pm \^ p.2 



(B8) 



3rd-order solutions: 



5(3)(fe;i) = 



H^k^d^k-^rl^ko r 

^^^^p ^fc(A;-fei23)^o(fei)^o(fe2)5o(fc3) [Dl{t)F^%{kuk2M) + C^{t)C,y^{kuk2M) 

+F2{t)Fsyyn{kl, k2, k3) + I3{t)lsym{kl,k2,k3) + J3{t)Jsym{ki, k2, ks) + /STs (t)i^sym (fcl , ^2 , ^s) 



e^^\k;t) = 



+L3{t)Lsym{kl,k2,k3) 

d^kid^k2 



(B9) 



Suik - ki23) 5o{ki)6o{k2)6o{k3) -{DiDj/H) Gi%{ki,k2, fes) - {C3/H) C,y^{ki,k2, fes) 



(27r)6 "-1^0/ "uv"'i/"uv"'z/"uv"o/ 1^ "sym 

-{{F3 - DiF2)/H}Fsymikuk2,k3) - {{h- D^F2)/H}h{t)hy^,(kiMM) 
-{Jzl H)Jsyui{ki,k2,k3) - {k3/H)Ksy^{ki,k2,k3) - (L3/ff)Lsym(fei, ^2, ^3) • 



(BIO) 
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The symmetrized kernels of integrals are 



GgL(fcl,fe2,fc3) - - 



2 ( 3 / 

— /3(fei,fc23) i fi{k2,kz) + -\a{k2,k3,) + a(fc3,fc2) 

1 f 5 / 

+— fe23) \ /3(fc2, fes) + 2 ("C^a, ^3) + a(A;3, fca) 

+ ^a(fc23,fci) |/3(fc2,fc3) + ^(a(A;2,fc3) + a{kz,k2) 
2 f 3 / 

— /3(fci,fc23) S P{k2,k3) + -[a{k2,k3) + a{k3,k2) 



Csym(fcl, ^2, ^3) 
Fsyin(fel,fc2,fe3) 
^ym(fel, fc2, fe3) 
Jsym{kl,k2,k3) 

Ksyniiki,k2,k3) 
isym(fel, fc2, ^3) 



1 f 5 / 

+ ^23) S /3(A;2, fc3) + 2 ("(^2, ^3) + a{k3,k2) 

+ -^a{k23,ki) |/3(A;2,fc3) + 1(^(^2, ^3) +a(A;3,A;2) 
[/3(fei, fe23) (1 - ^2,3) + (cyclic perm.)] , 
[Q;(fei, fe23) (1 - M2,3) + (cyclic perm.)] , 
[a(fe23, fci) (1 - /i2,3) + (cyclic perm.)] , 
[(1 - lA.2'i)Pik2,k3) + (cyclic perm.)] , 
^(1 - Aii,23){"('^2, ^3) + a(fc3, fe3)| + (cyclic perm 
[(1 - Aii,23)(l - A^L) + (cyclic perm.)] , 



(cyclic perm.) 



(cyclic perm.) 



In the expressions (|B9p and (jB10[) . new growth functions originating from non-linearity of the Poisson equation 
appear. The evolution equations for these function are given by 

CC3 = D1F2, 

Cl3^ Di{F2 + 2HF2)+DiF2, 



CJ3 -- 
CK3 



6/33 V 3 y 7 

6/33 V 3 J 7 ' 



2. One-loop power spectra 



Using the perturbative solutions obtained above, the power spectrum can be expressed as 
Pah{k;t)^Pll'\k;t)+P^l'\k;t)+P^^l'\k;t) + --- ; (a, 6 = ,5 or 
The terms P^^^\k) imply the linear power spectrum, given by 



P^l'\k;t) = Dlit)Po{k), pai)(fc;i) = -MlMp„(fc), P^l'\k;t)= £0-] Poik), 



H 

where the power spectrum Po(fc) is defined by 

(5o(fe)'5o(fc')) = (27r)3jD(fc + k') Po(|fc|). 



H 



(Bll) 



(B12) 



(B13) 
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The terms P^'^'^\k) and P^^^\k) are the so-called one-loop power spectra whose explicit expressions respectively 
become 



x2 



^1 1 14x(l + a;2 - 2mx) j + ^ ^ ^2 _ 2^3; 



(3X + 7/X- 10/i2x)(l -^2) 



7x{l + a;2 - 2/ix)2 



fc3 
(27r)^ 



DiDf {x-7iJ. + 6fi^x){-3x-7fi + WlJ.'^x) F2F2 2{ij,'^ - 1^ 



H 98x2(1 + a;2 - 2/i.T)2 H {1 + x"^ - 2fix)'^ 

F2DI {-2,x - 7/i + 10Ai2x)(/i2 _ 1) £)i£)iF2 (a; - T^i + 6Ai2a;)(^2 _ 



H 

(27r)- ^0 
x2 



7a;(l + x2 - 2/xa;)2 7a;(l + a;2 - 2/ia;)2 

dxx'^Po{kx) J d/LtPo ^fc-\/r+"^2rr2^j 



+ 



f {x-7f, + 6fj,''x) y F2\ (/x' - 1)' 
H j \Ux{l + x'^ -2nx)j \H J (1 + a;2 - 2/xa;)2 

D1F2D1 {x-7n + Qi?x){ii'^ - 1) 



if2 



7a;(l + x2 - 2nxf 



for P(22)(fc), and 



3 /■CO 

^Poik) / dxPoikx) 

Df{t) ^^^{l2a; - 158a;^ + lOOa;^ - 42a;^ + 3{x'^ - lf{7x^ + 2) In | 



+D^{t)Cz{t) -X;[Qx- l&x^ - Qx^ + 'Mx^ - if In 
+£)i(i)/3(i) ^{Sa; + 16a;^ - 6a;^ -h 3(a;2 - l)^ In 



l + x 



l-x 



} 



l-x 



} 



(B14) 



(B15) 



(B16) 



+D^{t){Ki{t)+Lr,{t)} ^{-6x + 22x3 + 22x5 -6x^ + 3(x2- 1)4 In }] , (B17) 



^3 /-oo 
JTTT^Poik) / dxPo(fcx) 

;tf^|24x - 202x3 + SSx^ - 30x^ + Z{x^ - 1)3(5x2 + 4) In \ 



H 



!l 1 |6x - 16x3 - ^ 3(^2 _ ^^3 1^ 1+^ 1 
12x'^ L 1 — X J 



H ^ I^{6^ + 16^^-6.- + 3(x2-l)3ln^} 

{Z^U/Vi + I,)}- 1 |_6a; + 22x3 + 22x5 -6x^ + 3(x2- 1)4 In 



24x3 



l + x 


}] 


l-x 





(BIS) 
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Pee "* 



^3 

J^Mk) I dxPoikx) 



DiDi 



^(Ux- 82x^ + Ax^ - + 3(x^ - lf(x^ + 2) In | 
Mx^ I 1 - X ) 



H 

H'^ ox'^ I 1 — X ) 

Di{h~F2Di) 1 



i/2 12x3 



■ (Gx + 16a;3 - + Sfa;^ - 1)^ In —- \ 

: L 1 — X ) 

|-6a; + 22a;3 + 22x^ - Qx'^ + 3(a;^ - 1)"^ In 



1 + a; 




1-a; 


}] 



(B19) 



for P^^^\k). 

The power spectrum without the non-Unear interaction terms X can be obtained by putting F2 — — — K3 
L3 = 0in P(22) and P(i3). 
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